﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SmartMathLibrary.LinearEquationSolvers;

namespace SmartMathLibrary.NonlinearEquationSolvers
{
    /// <summary>
    /// This class provides the solving of a nonlinear system of equation by using the extended
    /// Newton method.
    /// </summary>
    [Serializable]
    public class NewtonNonLinearSystemSolver : AbstractRealNonLinearSystemSolver
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="NewtonNonLinearSystemSolver"/> class.
        /// </summary>
        /// <param name="equationSystem">The nonlinear system of equation to solve.</param>
        public NewtonNonLinearSystemSolver(NonLinearEquationSystem equationSystem)
            : base(equationSystem)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="NewtonNonLinearSystemSolver"/> class.
        /// </summary>
        /// <param name="equationSystemFunctions">The nonlinear system of equation to solve.</param>
        public NewtonNonLinearSystemSolver(HardMultivariateRealFunction[] equationSystemFunctions)
            : base(equationSystemFunctions)
        {
        }

        /// <summary>
        /// Solves the nonlinear system of equation by using the extended Newton method. 
        /// The property PrecisionError have to be checked after the calculation to be sure 
        /// that the results are correct.
        /// </summary>
        /// <param name="startValues">The start values for the calculation.</param>
        /// <param name="iterations">The maximum number of iterations.</param>
        /// <param name="precision">The precision of the result.</param>
        /// <returns>
        /// The approximated numeric results, which solves the nonlinear system of equation.
        /// </returns>
        public override GeneralVector SolveNonLinearSystem(GeneralVector startValues, int iterations, double precision)
        {
            if (startValues == (GeneralVector) null)
            {
                throw new ArgumentNullException("startValues");
            }

            if (startValues.Count != this.EquationSystem.EquationSystem.Length)
            {
                throw new ArgumentException(
                    "The number of components in the start vector has to be even with the number of equations.");
            }

            if (this.EquationSystem.EquationSystem.Length < 2)
            {
                throw new ArgumentException(
                    "The system of equation has to have to contain a minimum number of two equations.");
            }

            int iterationCount = 0;
            const double epsilon = 1e-15;
            GeneralVector result = startValues.Copy();
            GeneralVector parameterVector = new GeneralVector(this.EquationSystem.EquationSystem.Length);
            Matrix jacobianMatrix = new Matrix(this.EquationSystem.EquationSystem.Length,
                                               this.EquationSystem.EquationSystem.Length);

            while ((iterationCount < iterations) && (!CheckPrecision(this.EquationSystem
                                                                         .SolveAt(result), precision)))
            {
                GeneralVector functionResults = this.EquationSystem.SolveAt(result);

                for (int i = 0; i < this.EquationSystem.EquationSystem.Length; i++)
                {
                    double delta = 0;
                    double tempuri = result[i];

                    if (tempuri > 1.0)
                    {
                        delta = epsilon * tempuri;
                    }
                    else
                    {
                        delta = epsilon;
                    }

                    result[i] = tempuri + delta;
                    delta = result[i] - tempuri;
                    parameterVector = this.EquationSystem.SolveAt(result);
                    result[i] = tempuri;

                    for (int j = 0; j < this.EquationSystem.EquationSystem.Length; j++)
                    {
                        jacobianMatrix.SetValueAtPosition(j, i, (parameterVector[j] - functionResults[j]) / delta);
                    }
                }

                for (int i = 0; i < this.EquationSystem.EquationSystem.Length; i++)
                {
                    double tempuri = 0.0;

                    for (int j = 0; j < this.EquationSystem.EquationSystem.Length; j++)
                    {
                        tempuri += jacobianMatrix.GetValueAtPosition(i, j) * result[j];
                    }

                    parameterVector[i] = tempuri - functionResults[i];
                }

                GaussLinearSystemSolver linearSystemSolver = new GaussLinearSystemSolver(jacobianMatrix, parameterVector);

                result = linearSystemSolver.SolveLinearSystem();
                iterationCount++;
            }

            this.neededIterations = iterationCount;
            this.precisionError = (iterations == iterationCount);

            return result;
        }

        /// <summary>
        /// Checks the precision of the calculation.
        /// </summary>
        /// <param name="checkVector">The vector to check.</param>
        /// <param name="precision">The precision to check.</param>
        /// <returns>True, if all components are smaller than the specified precision otherwise, 
        /// false.</returns>
        private static bool CheckPrecision(GeneralVector checkVector, double precision)
        {
            for (int i = 0; i < checkVector.Count; i++)
            {
                if (Math.Abs(checkVector[i]) > precision)
                {
                    return false;
                }
            }

            return true;
        }
    }
}